HerbPt1 16S Figures

require(phyloseq)
require(tidyverse)
require(phyloseq)
require(reshape2)

Loading required package: reshape2

Attaching package: ‘reshape2’

The following objects are masked from ‘package:data.table’:

    dcast, melt

The following object is masked from ‘package:tidyr’:

    smiths

require(dplyr)
require(ggplot2)

Load data

ps_dmn <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/DMN_ests_16S.Rdata")
sample_data(ps_dmn)$Herbicide <- factor(sample_data(ps_dmn)$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
ps_rare <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/HerbPt1_rare_16S.Rdata")
sample_data(ps_rare)$Herbicide <- factor(sample_data(ps_rare)$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
ps_trans <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/HerbPt1_hel_trans_16S.Rdata")
sample_data(ps_trans)$Herbicide <- factor(sample_data(ps_trans)$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

create alphadiversity tables

alpha_div <- estimate_richness(physeq = ps_rare, measures = c("Observed", "Shannon", "Chao1"))
#pull out metadata and concatonate with alpha diversity metrics
md<-data.frame(sample_data(ps_rare))
alpha_div_md <- rownames_to_column(alpha_div, "Barcode_ID_G") %>% full_join(md) 

Joining, by = "Barcode_ID_G"

alpha_div_md$Herbicide <- factor(alpha_div_md$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

Shannon Div plots - no significant differences among herbicide treatments at any of the three time points

ggplot(data = alpha_div_md, aes(Herbicide, Shannon, color= Herbicide)) + facet_grid(. ~ Time) + geom_boxplot() + theme_classic() + theme(axis.text.x = element_text(angle = 45, hjust = 1) )

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_Shannon.pdf")

Saving 7.29 x 4.51 in image

aov_t1<-aov(Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T1",])
plot(aov_t1$residuals)

summary(aov_t1)

            Df Sum Sq Mean Sq F value Pr(>F)
Herbicide    4 0.0689 0.01723   1.345  0.266
Residuals   51 0.6533 0.01281               

aov_t2<-aov(Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T2",])
plot(aov_t2$residuals)

summary(aov_t2)

            Df Sum Sq Mean Sq F value Pr(>F)
Herbicide    4 0.0449 0.01123   0.865  0.491
Residuals   49 0.6355 0.01297               

aov_t3<-aov(Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T3",])
plot(aov_t3$residuals)

summary(aov_t3)

            Df Sum Sq Mean Sq F value Pr(>F)
Herbicide    4 0.0457 0.01142   0.681  0.609
Residuals   50 0.8391 0.01678               

ordinations and adonis testing with three separate objects (i.e., dmn, rarefied, transformed). Rare taxa are removed from rarefied and transfomred to sucessfully ordinate. At this point, the transformed data will not ordinate.

ord_dmn<-ordinate(physeq = ps_dmn, method = "NMDS", distance = "bray", k=3, trymax= 100)

ps_rare_sub<-prune_taxa(taxa_sums(ps_rare) > 2, ps_rare)
ord_rare<-ordinate(physeq = ps_rare_sub, method = "NMDS", distance = "bray", k=3, trymax= 100, previous.best = ord_rare)

#can't get the hellinger transformed data to ordinate successfully
ps_trans_sub<-prune_taxa(taxa_sums(ps_trans) > 0.5, ps_trans)
ord_transformed<-ordinate(physeq = ps_trans_sub, method = "NMDS", distance = "bray", k=3, trymax= 100)

Adonis testing of herbicide treatments by time point

ps_adonis<-function(physeq){
  otu_tab<-data.frame(phyloseq::otu_table(physeq))
  md_tab<-data.frame(phyloseq::sample_data(physeq))
    if(taxa_are_rows(physeq)== T){
       physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
            else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
  print(anova(vegan::betadisper(physeq_dist, md_tab$Herbicide)))
  vegan::adonis(physeq_dist ~ Herbicide * Time, data = md_tab)
}

ps_adonis(ps_rare_sub)

Analysis of Variance Table

Response: Distances
           Df   Sum Sq    Mean Sq F value Pr(>F)
Groups      4 0.002281 0.00057033  0.4454 0.7756
Residuals 160 0.204866 0.00128041               

Call:
vegan::adonis(formula = physeq_dist ~ Herbicide * Time, data = md_tab) 

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

                Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
Herbicide        4    0.7552 0.18880  1.1017 0.02657  0.037 * 
Time             2    0.3792 0.18959  1.1063 0.01334  0.060 . 
Herbicide:Time   8    1.5797 0.19747  1.1523 0.05559  0.003 **
Residuals      150   25.7059 0.17137         0.90450          
Total          164   28.4200                 1.00000          
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

ps_adonis(ps_trans_sub)

Analysis of Variance Table

Response: Distances
           Df  Sum Sq   Mean Sq F value Pr(>F)
Groups      4 0.01133 0.0028322  0.3169 0.8664
Residuals 165 1.47460 0.0089369               

Call:
vegan::adonis(formula = physeq_dist ~ Herbicide * Time, data = md_tab) 

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

                Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)  
Herbicide        4    0.2416 0.060388 0.97118 0.02277  0.573  
Time             2    0.1771 0.088564 1.42433 0.01670  0.012 *
Herbicide:Time   8    0.5526 0.069078 1.11094 0.05209  0.067 .
Residuals      155    9.6378 0.062179         0.90845         
Total          169   10.6091                  1.00000         
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

ps_adonis(ps_dmn)

Analysis of Variance Table

Response: Distances
           Df   Sum Sq   Mean Sq F value  Pr(>F)  
Groups      4 0.019549 0.0048873  2.8407 0.02596 *
Residuals 165 0.283879 0.0017205                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Call:
vegan::adonis(formula = physeq_dist ~ Herbicide * Time, data = md_tab) 

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

                Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)    
Herbicide        4    0.2226 0.055638  2.2900 0.04823  0.001 ***
Time             2    0.1044 0.052199  2.1484 0.02263  0.008 ** 
Herbicide:Time   8    0.5214 0.065169  2.6823 0.11299  0.001 ***
Residuals      155    3.7659 0.024296         0.81615           
Total          169    4.6142                  1.00000           
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ordination plots DMN

ord_t1_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Run 0 stress 0.0860376 
Run 1 stress 0.08603687 
... New best solution
... Procrustes: rmse 0.0006267826  max resid 0.003492555 
... Similar to previous best
Run 2 stress 0.08603727 
... Procrustes: rmse 0.000519486  max resid 0.002879424 
... Similar to previous best
Run 3 stress 0.08603775 
... Procrustes: rmse 0.0006710548  max resid 0.00365014 
... Similar to previous best
Run 4 stress 0.08674028 
Run 5 stress 0.08674051 
Run 6 stress 0.08610309 
... Procrustes: rmse 0.0155372  max resid 0.1010944 
Run 7 stress 0.08662927 
Run 8 stress 0.08674034 
Run 9 stress 0.08603713 
... Procrustes: rmse 0.0004570747  max resid 0.002500815 
... Similar to previous best
Run 10 stress 0.1023472 
Run 11 stress 0.08674014 
Run 12 stress 0.09941847 
Run 13 stress 0.08610274 
... Procrustes: rmse 0.01547852  max resid 0.1009111 
Run 14 stress 0.1017673 
Run 15 stress 0.08662931 
Run 16 stress 0.08610275 
... Procrustes: rmse 0.01547835  max resid 0.1009027 
Run 17 stress 0.08662977 
Run 18 stress 0.08662964 
Run 19 stress 0.08603738 
... Procrustes: rmse 0.0005091999  max resid 0.002807653 
... Similar to previous best
Run 20 stress 0.08610325 
... Procrustes: rmse 0.01555667  max resid 0.1011618 
*** Solution reached

T1_dmn<-ggordiplots::gg_ordiplot(ord = ord_t1_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T1_dmn$plot + theme_classic()

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T1.pdf")

Saving 7.29 x 4.51 in image

ord_t2_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T2"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Run 0 stress 0.1104159 
Run 1 stress 0.1104149 
... New best solution
... Procrustes: rmse 0.001171513  max resid 0.007161514 
... Similar to previous best
Run 2 stress 0.1104149 
... Procrustes: rmse 0.0002875201  max resid 0.001118723 
... Similar to previous best
Run 3 stress 0.109526 
... New best solution
... Procrustes: rmse 0.03334389  max resid 0.2369979 
Run 4 stress 0.1095359 
... Procrustes: rmse 0.0009442239  max resid 0.004563832 
... Similar to previous best
Run 5 stress 0.1121088 
Run 6 stress 0.1095264 
... Procrustes: rmse 0.0001397644  max resid 0.0007317152 
... Similar to previous best
Run 7 stress 0.1122111 
Run 8 stress 0.1095263 
... Procrustes: rmse 0.0008604593  max resid 0.003402465 
... Similar to previous best
Run 9 stress 0.1095262 
... Procrustes: rmse 0.0008348408  max resid 0.003298233 
... Similar to previous best
Run 10 stress 0.1095262 
... Procrustes: rmse 0.0005129767  max resid 0.00314258 
... Similar to previous best
Run 11 stress 0.1108775 
Run 12 stress 0.1247908 
Run 13 stress 0.1104148 
Run 14 stress 0.1095266 
... Procrustes: rmse 0.0009426628  max resid 0.004072489 
... Similar to previous best
Run 15 stress 0.1122114 
Run 16 stress 0.1122109 
Run 17 stress 0.1104146 
Run 18 stress 0.1120888 
Run 19 stress 0.1120856 
Run 20 stress 0.1108786 
*** Solution reached

T2_dmn<-ggordiplots::gg_ordiplot(ord = ord_t2_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T2_dmn$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T2.pdf")

Saving 7.29 x 4.51 in image

ord_t3_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T3"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Run 0 stress 0.09696447 
Run 1 stress 0.09982855 
Run 2 stress 0.09784799 
Run 3 stress 0.09800983 
Run 4 stress 0.09784894 
Run 5 stress 0.0969634 
... New best solution
... Procrustes: rmse 0.04061371  max resid 0.2635214 
Run 6 stress 0.09706969 
... Procrustes: rmse 0.04032592  max resid 0.2678095 
Run 7 stress 0.09696485 
... Procrustes: rmse 0.04056564  max resid 0.2625302 
Run 8 stress 0.0969324 
... New best solution
... Procrustes: rmse 0.004791775  max resid 0.03220976 
Run 9 stress 0.09696451 
... Procrustes: rmse 0.03977923  max resid 0.263354 
Run 10 stress 0.09696596 
... Procrustes: rmse 0.004560318  max resid 0.02412528 
Run 11 stress 0.09784882 
Run 12 stress 0.09838599 
Run 13 stress 0.1005549 
Run 14 stress 0.09777523 
Run 15 stress 0.0969219 
... New best solution
... Procrustes: rmse 0.001739056  max resid 0.008015476 
... Similar to previous best
Run 16 stress 0.1008367 
Run 17 stress 0.09697109 
... Procrustes: rmse 0.03948731  max resid 0.2623541 
Run 18 stress 0.09696141 
... Procrustes: rmse 0.005146037  max resid 0.03412768 
Run 19 stress 0.1028377 
Run 20 stress 0.09707004 
... Procrustes: rmse 0.03882214  max resid 0.2583721 
*** Solution reached

T3_dmn<-ggordiplots::gg_ordiplot(ord = ord_t3_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T3_dmn$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T3.pdf")

Saving 7.29 x 4.51 in image

Ordination plots rarefied

ord_t1_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1810491 
Run 1 stress 0.181593 
Run 2 stress 0.191872 
Run 3 stress 0.1815709 
Run 4 stress 0.1820694 
Run 5 stress 0.1815727 
Run 6 stress 0.1814585 
... Procrustes: rmse 0.0482874  max resid 0.155698 
Run 7 stress 0.1831936 
Run 8 stress 0.1821565 
Run 9 stress 0.1811472 
... Procrustes: rmse 0.02328168  max resid 0.1241755 
Run 10 stress 0.1809422 
... New best solution
... Procrustes: rmse 0.01257791  max resid 0.07384735 
Run 11 stress 0.1830541 
Run 12 stress 0.1821389 
Run 13 stress 0.1819163 
Run 14 stress 0.1813678 
... Procrustes: rmse 0.01983823  max resid 0.1171532 
Run 15 stress 0.1821569 
Run 16 stress 0.1830532 
Run 17 stress 0.1819142 
Run 18 stress 0.1837675 
Run 19 stress 0.1860978 
Run 20 stress 0.1811476 
... Procrustes: rmse 0.01835014  max resid 0.1232744 
Run 21 stress 0.1826836 
Run 22 stress 0.182507 
Run 23 stress 0.1817998 
Run 24 stress 0.1825066 
Run 25 stress 0.1815422 
Run 26 stress 0.1819123 
Run 27 stress 0.1816038 
Run 28 stress 0.181148 
... Procrustes: rmse 0.0183826  max resid 0.1233536 
Run 29 stress 0.1915283 
Run 30 stress 0.1823893 
Run 31 stress 0.1821972 
Run 32 stress 0.1827478 
Run 33 stress 0.1837639 
Run 34 stress 0.1823302 
Run 35 stress 0.1856634 
Run 36 stress 0.1837639 
Run 37 stress 0.1811488 
... Procrustes: rmse 0.01895138  max resid 0.1248699 
Run 38 stress 0.180967 
... Procrustes: rmse 0.02557771  max resid 0.1286128 
Run 39 stress 0.183194 
Run 40 stress 0.1837639 
Run 41 stress 0.181913 
Run 42 stress 0.1825652 
Run 43 stress 0.1842313 
Run 44 stress 0.1825074 
Run 45 stress 0.1819123 
Run 46 stress 0.1830533 
Run 47 stress 0.1810718 
... Procrustes: rmse 0.02863214  max resid 0.1313275 
Run 48 stress 0.1917654 
Run 49 stress 0.1831936 
Run 50 stress 0.1834936 
Run 51 stress 0.192229 
Run 52 stress 0.182276 
Run 53 stress 0.1822762 
Run 54 stress 0.1825078 
Run 55 stress 0.1815184 
Run 56 stress 0.1821958 
Run 57 stress 0.1814159 
... Procrustes: rmse 0.04113413  max resid 0.1452096 
Run 58 stress 0.1814153 
... Procrustes: rmse 0.04101125  max resid 0.145104 
Run 59 stress 0.186166 
Run 60 stress 0.1814141 
... Procrustes: rmse 0.0407537  max resid 0.1447673 
Run 61 stress 0.181331 
... Procrustes: rmse 0.03588519  max resid 0.1397534 
Run 62 stress 0.1822761 
Run 63 stress 0.1808178 
... New best solution
... Procrustes: rmse 0.01546488  max resid 0.09409532 
Run 64 stress 0.1811484 
... Procrustes: rmse 0.02374035  max resid 0.1223041 
Run 65 stress 0.1822021 
Run 66 stress 0.1822765 
Run 67 stress 0.1814159 
Run 68 stress 0.1813262 
Run 69 stress 0.1825039 
Run 70 stress 0.1816037 
Run 71 stress 0.1816037 
Run 72 stress 0.1823672 
Run 73 stress 0.182367 
Run 74 stress 0.1810519 
... Procrustes: rmse 0.02345899  max resid 0.1072637 
Run 75 stress 0.181525 
Run 76 stress 0.1825929 
Run 77 stress 0.1810721 
... Procrustes: rmse 0.03648172  max resid 0.1453981 
Run 78 stress 0.1813256 
Run 79 stress 0.1880286 
Run 80 stress 0.1874007 
Run 81 stress 0.185664 
Run 82 stress 0.1814029 
Run 83 stress 0.1808177 
... New best solution
... Procrustes: rmse 0.000342235  max resid 0.001823038 
... Similar to previous best
*** Solution reached

T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T1_rare$plot + theme_classic()

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T1.pdf")

Saving 7.29 x 4.51 in image

ord_t2_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T2"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1534021 
Run 1 stress 0.1512024 
... New best solution
... Procrustes: rmse 0.0754476  max resid 0.2676264 
Run 2 stress 0.1493358 
... New best solution
... Procrustes: rmse 0.06201399  max resid 0.2214984 
Run 3 stress 0.1532387 
Run 4 stress 0.1493338 
... New best solution
... Procrustes: rmse 0.03039888  max resid 0.1722904 
Run 5 stress 0.1505376 
Run 6 stress 0.1508442 
Run 7 stress 0.1538361 
Run 8 stress 0.1493389 
... Procrustes: rmse 0.02035211  max resid 0.1193448 
Run 9 stress 0.1493331 
... New best solution
... Procrustes: rmse 0.03006402  max resid 0.1666307 
Run 10 stress 0.1527413 
Run 11 stress 0.1494084 
... Procrustes: rmse 0.02890853  max resid 0.1471386 
Run 12 stress 0.1503514 
Run 13 stress 0.1564968 
Run 14 stress 0.1493168 
... New best solution
... Procrustes: rmse 0.01404071  max resid 0.07782122 
Run 15 stress 0.1493021 
... New best solution
... Procrustes: rmse 0.01741756  max resid 0.09981203 
Run 16 stress 0.1493198 
... Procrustes: rmse 0.02215659  max resid 0.1213156 
Run 17 stress 0.1507341 
Run 18 stress 0.149334 
... Procrustes: rmse 0.00663877  max resid 0.03659659 
Run 19 stress 0.1498004 
... Procrustes: rmse 0.03768561  max resid 0.1683317 
Run 20 stress 0.1495337 
... Procrustes: rmse 0.03069104  max resid 0.1677622 
Run 21 stress 0.1505776 
Run 22 stress 0.1493306 
... Procrustes: rmse 0.02690526  max resid 0.1437716 
Run 23 stress 0.1516113 
Run 24 stress 0.1498424 
Run 25 stress 0.1493193 
... Procrustes: rmse 0.02155654  max resid 0.1168637 
Run 26 stress 0.1508459 
Run 27 stress 0.1497011 
... Procrustes: rmse 0.03123585  max resid 0.1580844 
Run 28 stress 0.1499895 
Run 29 stress 0.1494448 
... Procrustes: rmse 0.007164553  max resid 0.03016139 
Run 30 stress 0.1498771 
Run 31 stress 0.1502641 
Run 32 stress 0.1558884 
Run 33 stress 0.1501174 
Run 34 stress 0.1498656 
Run 35 stress 0.1508896 
Run 36 stress 0.1520205 
Run 37 stress 0.1495492 
... Procrustes: rmse 0.01320733  max resid 0.04627731 
Run 38 stress 0.1493166 
... Procrustes: rmse 0.002633167  max resid 0.01254371 
Run 39 stress 0.1499219 
Run 40 stress 0.1554524 
Run 41 stress 0.1508795 
Run 42 stress 0.1529847 
Run 43 stress 0.1501535 
Run 44 stress 0.1496195 
... Procrustes: rmse 0.04649384  max resid 0.2115397 
Run 45 stress 0.1502646 
Run 46 stress 0.1518985 
Run 47 stress 0.1497692 
... Procrustes: rmse 0.03326784  max resid 0.1546212 
Run 48 stress 0.1509269 
Run 49 stress 0.1502629 
Run 50 stress 0.1533299 
Run 51 stress 0.1544682 
Run 52 stress 0.1541878 
Run 53 stress 0.1541584 
Run 54 stress 0.1501233 
Run 55 stress 0.1501996 
Run 56 stress 0.1492999 
... New best solution
... Procrustes: rmse 0.001417234  max resid 0.006130985 
... Similar to previous best
*** Solution reached

T2_rare<-ggordiplots::gg_ordiplot(ord = ord_t2_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T2_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T2.pdf")

Saving 7.29 x 4.51 in image

ord_t3_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T3"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1249587 
Run 1 stress 0.1256869 
Run 2 stress 0.1257678 
Run 3 stress 0.127937 
Run 4 stress 0.1278503 
Run 5 stress 0.124948 
... New best solution
... Procrustes: rmse 0.00253675  max resid 0.00957671 
... Similar to previous best
Run 6 stress 0.1258153 
Run 7 stress 0.1256208 
Run 8 stress 0.1257204 
Run 9 stress 0.1255412 
Run 10 stress 0.1250279 
... Procrustes: rmse 0.006195166  max resid 0.03810374 
Run 11 stress 0.1250437 
... Procrustes: rmse 0.007767753  max resid 0.03872528 
Run 12 stress 0.1258564 
Run 13 stress 0.1249904 
... Procrustes: rmse 0.01078475  max resid 0.03446339 
Run 14 stress 0.1256788 
Run 15 stress 0.1251965 
... Procrustes: rmse 0.01742177  max resid 0.0711806 
Run 16 stress 0.1258553 
Run 17 stress 0.1249628 
... Procrustes: rmse 0.009053109  max resid 0.02873361 
Run 18 stress 0.1256671 
Run 19 stress 0.1249994 
... Procrustes: rmse 0.0125122  max resid 0.04112397 
Run 20 stress 0.1257192 
*** Solution reached

T3_rare<-ggordiplots::gg_ordiplot(ord = ord_t3_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T3_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T3.pdf")

Saving 7.29 x 4.51 in image

box and whisker plots of distance within group distances

library(micrUBIfuns)
beta_boxplot(physeq = subset_samples(ps_rare, Time=="T1"), method = "bray", group = "Herbicide")

$data

$plot

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T1_rare_withingroup_beta.pdf")

Saving 7.29 x 4.51 in image

beta_boxplot(physeq = subset_samples(ps_rare, Time=="T2"), method = "bray", group = "Herbicide")

$data

$plot

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T2_rare_withingroup_beta.pdf")

Saving 7.29 x 4.51 in image

beta_boxplot(physeq = subset_samples(ps_rare, Time=="T3"), method = "bray", group = "Herbicide")

$data

$plot

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T3_rare_withingroup_beta.pdf")

Saving 7.29 x 4.51 in image

treatment to control

plotDistances = function(p, m, s, d) {

  # calc distances
  wu = phyloseq::distance(p, m)
  wu.m = melt(as.matrix(wu))
  
  # remove self-comparisons
  wu.m = wu.m %>%
    filter(as.character(Var1) != as.character(Var2)) %>%
    mutate_if(is.factor,as.character)
  
  # get sample data (S4 error OK and expected)
  sd = data.frame(sample_data(p)) %>%
    select(s, d) %>%
    mutate_if(is.factor,as.character)
  sd$Herbicide <- factor(sd$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
  
  # combined distances with sample data
  colnames(sd) = c("Var1", "Type1")
  wu.sd = left_join(wu.m, sd, by = "Var1")
  
  colnames(sd) = c("Var2", "Type2")
  wu.sd = left_join(wu.sd, sd, by = "Var2")
  
  #remove this line to plot all comparisons. 
  wu.sd = wu.sd %>% filter(Type1 == "Hand" | Type1 == "Non-Treated")
  
  # plot
  ggplot(wu.sd, aes(x = Type2, y = value)) +
    theme_bw() +
    geom_point() +
    geom_boxplot(aes(color = ifelse(Type1 == Type2, "red", "black"))) +
    scale_color_identity() +
    facet_wrap(~ Type1, scales = "free_x") +
    theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
    ggtitle(paste0("Distance Metric = ", m))
  
}

a<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T1"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
a <- a + ggtitle("Time 1 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T1_rare_allgroup_beta.pdf")

Saving 7 x 7 in image

b<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T2"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
b <-b + ggtitle("Time 2 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T2_rare_allgroup_beta.pdf")

Saving 7 x 7 in image

c<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T3"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
c<- c + ggtitle("Time 3 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T3_rare_allgroup_beta.pdf")

Saving 7 x 7 in image

library(ggpubr)
ggarrange(a, b, c, ncol = 1)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_allgroup_beta.pdf", width = 7, height = 10)

Taxon abundance bar plot

#create super long color vector
col_vector <- c("#000000", "#FFFF00", "#1CE6FF", "#FF34FF", "#FF4A46", "#008941", "#006FA6", "#A30059",
        "#FFDBE5", "#7A4900", "#0000A6", "#63FFAC", "#B79762", "#004D43", "#8FB0FF", "#997D87",
        "#5A0007", "#809693", "#FEFFE6", "#1B4400", "#4FC601", "#3B5DFF", "#4A3B53", "#FF2F80",
        "#61615A", "#BA0900", "#6B7900", "#00C2A0", "#FFAA92", "#FF90C9", "#B903AA", "#D16100",
        "#DDEFFF", "#000035", "#7B4F4B", "#A1C299", "#300018", "#0AA6D8", "#013349", "#00846F",
        "#372101", "#FFB500", "#C2FFED", "#A079BF", "#CC0744", "#C0B9B2", "#C2FF99", "#001E09",
        "#00489C", "#6F0062", "#0CBD66", "#EEC3FF", "#456D75", "#B77B68", "#7A87A1", "#788D66",
        "#885578", "#FAD09F", "#FF8A9A", "#D157A0", "#BEC459", "#456648", "#0086ED", "#886F4C",
        
        "#34362D", "#B4A8BD", "#00A6AA", "#452C2C", "#636375", "#A3C8C9", "#FF913F", "#938A81",
        "#575329", "#00FECF", "#B05B6F", "#8CD0FF", "#3B9700", "#04F757", "#C8A1A1", "#1E6E00",
        "#7900D7", "#A77500", "#6367A9", "#A05837", "#6B002C", "#772600", "#D790FF", "#9B9700",
        "#549E79", "#FFF69F", "#201625", "#72418F", "#BC23FF", "#99ADC0", "#3A2465", "#922329",
        "#5B4534", "#FDE8DC", "#404E55", "#0089A3", "#CB7E98", "#A4E804", "#324E72", "#6A3A4C",
        "#83AB58", "#001C1E", "#D1F7CE", "#004B28", "#C8D0F6", "#A3A489", "#806C66", "#222800",
        "#BF5650", "#E83000", "#66796D", "#DA007C", "#FF1A59", "#8ADBB4", "#1E0200", "#5B4E51",
        "#C895C5", "#320033", "#FF6832", "#66E1D3", "#CFCDAC", "#D0AC94", "#7ED379", "#012C58",
        
        "#7A7BFF", "#D68E01", "#353339", "#78AFA1", "#FEB2C6", "#75797C", "#837393", "#943A4D",
        "#B5F4FF", "#D2DCD5", "#9556BD", "#6A714A", "#001325", "#02525F", "#0AA3F7", "#E98176",
        "#DBD5DD", "#5EBCD1", "#3D4F44", "#7E6405", "#02684E", "#962B75", "#8D8546", "#9695C5",
        "#E773CE", "#D86A78", "#3E89BE", "#CA834E", "#518A87", "#5B113C", "#55813B", "#E704C4",
        "#00005F", "#A97399", "#4B8160", "#59738A", "#FF5DA7", "#F7C9BF", "#643127", "#513A01",
        "#6B94AA", "#51A058", "#A45B02", "#1D1702", "#E20027", "#E7AB63", "#4C6001", "#9C6966",
        "#64547B", "#97979E", "#006A66", "#391406", "#F4D749", "#0045D2", "#006C31", "#DDB6D0",
        "#7C6571", "#9FB2A4", "#00D891", "#15A08A", "#BC65E9", "#FFFFFE", "#C6DC99", "#203B3C",

        "#671190", "#6B3A64", "#F5E1FF", "#FFA0F2", "#CCAA35", "#374527", "#8BB400", "#797868",
        "#C6005A", "#3B000A", "#C86240", "#29607C", "#402334", "#7D5A44", "#CCB87C", "#B88183",
        "#AA5199", "#B5D6C3", "#A38469", "#9F94F0", "#A74571", "#B894A6", "#71BB8C", "#00B433",
        "#789EC9", "#6D80BA", "#953F00", "#5EFF03", "#E4FFFC", "#1BE177", "#BCB1E5", "#76912F",
        "#003109", "#0060CD", "#D20096", "#895563", "#29201D", "#5B3213", "#A76F42", "#89412E",
        "#1A3A2A", "#494B5A", "#A88C85", "#F4ABAA", "#A3F3AB", "#00C6C8", "#EA8B66", "#958A9F",
        "#BDC9D2", "#9FA064", "#BE4700", "#658188", "#83A485", "#453C23", "#47675D", "#3A3F00",
        "#061203", "#DFFB71", "#868E7E", "#98D058", "#6C8F7D", "#D7BFC2", "#3C3E6E", "#D83D66",

        "#2F5D9B", "#6C5E46", "#D25B88", "#5B656C", "#00B57F", "#545C46", "#866097", "#365D25",
        "#252F99", "#00CCFF", "#674E60", "#FC009C", "#92896B")

phylumGlommed <- tax_glom(ps_rare, "Phylum")

#t1
phylumGlommed_herb_t1 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T1"), group = "Herbicide")

Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion

phylumGlommed_herb_t1 <- transform_sample_counts(phylumGlommed_herb_t1, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t1)$Herbicide <- factor(sample_data(phylumGlommed_herb_t1)$Herbicide, levels = c(1, 2, 3, 4, 5), 
       labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

plot_bar(phylumGlommed_herb_t1, x = "Herbicide", fill = "Phylum")  + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/Taxon_barplot_t1.pdf")

Saving 7.29 x 4.51 in image

#t2
phylumGlommed_herb_t2 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T2"), group = "Herbicide")

Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion

phylumGlommed_herb_t2 <- transform_sample_counts(phylumGlommed_herb_t2, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t2)$Herbicide <- factor(sample_data(phylumGlommed_herb_t2)$Herbicide, levels = c(1, 2, 3, 4, 5), 
       labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

plot_bar(phylumGlommed_herb_t2, x = "Herbicide", fill = "Phylum")  + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/Taxon_barplot_t2.pdf")

Saving 7.29 x 4.51 in image

#t3
phylumGlommed_herb_t3 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T3"), group = "Herbicide")

Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion

phylumGlommed_herb_t3 <- transform_sample_counts(phylumGlommed_herb_t3, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t3)$Herbicide <- factor(sample_data(phylumGlommed_herb_t3)$Herbicide, levels = c(1, 2, 3, 4, 5), 
       labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

plot_bar(phylumGlommed_herb_t3, x = "Herbicide", fill = "Phylum")  + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/Taxon_barplot_t3.pdf")

Saving 7.29 x 4.51 in image